BINARY INTERACTION ALGORITHMS FOR THE SIMULATION OF 
FLOCKING AND SWARMING DYNAMICS* 

G. ALBlt AND L. PARESCfflt 

Abstract. Microscopic models of flocking and swarming takes in account large numbers of interacting individ- 
uals. Numerical resolution of large flocks implies huge computational costs. Typically for N interacting individuals 
we have a cost of 0(N 2 ). We tackle the problem numerically by considering approximated binary interaction dy- 

^_l namics described by kinetic equations and simulating such equations by suitable stochastic methods. This approach 

permits to compute approximate solutions as functions of a small scaling parameter e at a reduced complexity of 

^O O(N) operations. Several numerical results show the efficiency of the algorithms proposed. 
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1. Introduction. The study of mathematical models describing collective behavior and syn- 
chronized motion of animals, like bird flocks, fish schools and insect swarms, has attracted a lot 
of attention in recent years [TJ [TTJ [16j [13j [371 HH1 EH ESJ [42] . In biological systems such behaviors 
are observed in every level of the food chain, from the swarm intelligence of the zoo plankton, 
to bird flocking and fish schools, to mammals moving in formation [24} [29j [38]. Beside biology, 
emerging collective behaviors play a relevant role in several applications involving the dynamics of 
a large number of individuals/particles which range from computer science [41], physics [26J and 
engineering [32J to social sciences and economy [10J. We refer to [37J for a recent review of some 
of the mathematical topics and the applications involved. 

Naturally occurring synchronized motion has inspired several directions of research within the 
control community. A well-known application is related to formation flying missions and missions 
involving the coordinated control of several autonomous vehicles [31J. There are several current 
projects which are dealing with the formation flying and coordinated control of satellites, like the 
DARWIN project of the European Space Agency (ESA) with the goal of launching a space-based 
telescope aiding in the search for possible life-supporting planets, or the PRISMA project led by 
^ the Swedish Space Corporation (SSC) which will be the first real formation flying space mission 

launched [19J. 

In this manuscript we will focus on general models which are capable to reproduce flocking, 
swarming and other collective behaviors. Most of the classical models describing these phenomena 
are based on the simple definition of three interacting zones, the so called three-zones model [3, 30J. 

Let us briefly summarize the three-zones assumptions. We define three regions around each 
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individual: a short-range repulsion zone, an intermediate velocity alignment zone and a long-range 
^-H attraction zone (see Figure [TTT] ) . Each interaction between individuals is evaluated accordingly to 

the relative position in the model. 

• Repulsion zone: when individuals are too close each other they tend to move away from 
^^ that area. 

• Alignment zone: individuals try to identify the possible direction of the group and to align 
with it. 

• Attraction zone: when individuals are too far from the group they want to get closer. 
Typically different interaction models are taken in the different zones [TJ [22] or the modelling is 

focused on a specific zone, like the alignment/consensus dynamic [18, 35J. Of course the particular 
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Figure 1.1. Sketch of the three-zone model 

shape and size of the zones depends on the specific application considered. For example recent 
studied on birds flocks suggest that each bird modifies its position, relative to few individuals 
directly surrounding it, no matter how close or how far away those individuals are |5J. It is not 
clear however if this applies also to other kind of animals. 

Studying this kind of dynamics for large system of individuals implies a considerable effort in 
numerical simulations, microscopic models based on real data my take into account very large num- 
bers of interacting individuals (from several hundred thousands up to millions). Computationally 
the problems have the same structure of many classical problems in computational physics which 
require the evaluation of all pairwise long range interactions in a large ensembles of particles. The 
TV-body problem of gravitation (or electrostatics) is a classical example. Such problems involve 
the evaluation of summations of the type 

N 

S ? = X>itf(*i,a;), Vi. (1.1) 

i=i 

A direct evaluation of such sums at TV target points clearly requires 0(N 2 ) operations and algo- 
rithms which reduce the cost to 0(N a ) with 1 < a < 2, 0(N log N) etc. are referred to as fast 
summation methods. For uniform grid data the most famous of these is certainly the Fast Fourier 
Transform (FFT). In the case of general data most fast summation methods are approximate meth- 
ods based on analytical considerations, like the Fast Multipole method [26J, Wavelets Transform 
methods |6J and, more recently, dimension reduction using Compressive Sampling techniques [10J, 
or based on some Monte Carlo strategies at different levels |9J . Extensions of the above mentioned 
approaches to kinetic equations are discussed for example in jSH [33j HJ [55j [35] • 

From a mathematical modeling point of view, these problem have been developed extensively 
in the kinetic research community (see [5UJ[MJ[TII]) where the derivation of kinetic and an hydro- 
dynamic equations represent a first step towards the reduction of the computational complexity. 
Of course, passing from a microscopic description based on phase-space particles (xi(t) , Vi(t)) to 
a mesoscopic level where the object of study is a particle distribution function f(x, v, t) redefines 
the model in a new one where new methods of solution are required. 

In this paper we are going to follow this research path in two main directions: first we review 
the derivation of the different kinetic approximations from the original microscopic models and 
then we introduce and analyze several stochastic Monte Carlo methods to approximate the kinetic 
equations. Monte Carlo methods are the most well-known approach for the numerical solution 
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of the Boltzmann equation of rarefied gases in the short-range interactions, and many efficient 
algorithms have been presented [3 HI [9, 39J. On the other hand the literature on efficient Monte 
Carlo strategies for long-range interactions, and thus Landau-Fokker-Plank equations, is much less 
developed but of great interest in the field of plasma physics [HI [H] • 

Here, inspired by the techniques introduced in [8, 21J for plasma, we develop direct simulation 
Monte Carlo methods based on a binary collision dynamic described by the corresponding kinetic 
equation. The methods permit to approximate the microscopic dynamic at a cost directly propor- 
tional to the number of sample particles involved in the computation, thus avoiding the quadratic 
computational cost. The limiting behavior characterizing the mean- field interaction process of the 
particles system is recovered under a suitable asymptotic scaling of the binary collision process. 
In such a limit we show that the Monte Carlo methods here developed are in very good agree- 
ment with the direct evaluation of the original microscopic model but with a considerable gain of 
computational efficiency. 

The rest of the manuscript is organized as follows. In the first section we present some of 
the classical microscopic models for flocking and swarming. Generalization of the notion of visual 
cone [13J are also discussed. Since the interaction is non local, the derivation of the limiting 
mean-field kinetic equation is made through a Povzner-Boltzmann kinetic equation [40J in the 
anologous situation of the so-called grazing collision limit [12J. To solve the resulting Boltzmann- 
like mesoscopic partial differential equation we introduce different stochastic binary interaction 
algorithms and compare their computational efficiency and accuracy with a direct evaluation of 
the microscopic models and a stochastic approximation of the mean-field kinetic model. We show 
that the new approach permits to reduce the overall cost from 0(N 2 ) to O(N) operations. In 
particular we show that the choice e = At, where e is the small scaling parameter leading to 
the mean field kinetic model, originates binary interaction algorithms consistent with the limiting 
behavior of the particle system. Furthermore, in contrast with classical methods [8, 21J, the nature 
of the approximating equations is such that the resulting Monte Carlo algorithms are fully mesh 
less. In the last section of the paper we report several simulations in two and three space dimensions 
of different microscopic models solved by the binary Monte Carlo method in the above scaling. 

2. Microscopic models. In this section we review some well-known microscopic models of 
flocking and swarming (see [18l|22j[35] and the references therein). We are interested in the study 
of a dynamical system composed of N individuals with the following general structure 

X{ = Vi 

! N i = l,...,7V, 

Vi = S(vi) + j~^2[H(xi,Xj)(vj -v^ +A(xi,Xj) + R(xi,Xj)]il) a (xi,Xj,Vi) 

3 = 1 

(2.1) 
where {xi,Vi) lives in R 2d , d > 1, S(vi) describes a self-propelling term, H(xi,Xj) the alignment 
process, A(xi, Xj) the attraction dynamic and the term R(xi,Xj) the short-range repulsion. In (2.1 ) 
the multiplicative factor iJj a (xi,Xj,Vi) G [0, 1] takes into account the effects of space perception as 
a function of some vector of parameters a. 

2.1. Cucker and Smale model. Cucker and Smale model is a pure alignment model, no 
repulsion or attraction or other effects are taken in account, see [T51 [T7] and [12J. The classical 
model reads as follow 



1 n i = l,...,JV, (2.2) 

3 = 1 
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Figure 2.1. Profile of the function H in Cucker-Smale model 



where H(\xj — Xi\) is a function that measures the strength of the interaction between individuals 
i and j , and depends on the mutual distance, under the assumption that closer individuals have 
more influence than the far distance ones. 

A typical choice of function H is the following 



H(r) 



K 



(c 2 + r 2 )7' 



(2.3) 



where K^ > are positive parameters and 7 > 0. Under this assumptions it can be shown that 



well-posedness holds for the initial value problem of (2.2) and the solution is mass and momentum 
preserving, with compact support for position and velocity, see [27J. 

Moreover in p~8| [T2] it was established that the parameter 7 discriminate the behavior of the 
solution, in the following way 

Theorem 2.1. Let T(i) = \Y.i± 3 - M*) - x j( t )\ 2 and K(t) = \\vi(t) - Vj(t)\ 2 . If 7 < \ then 
(i) exist a positive costant Bq such that: T(t) < Bq for all t G R. 
(ii) A(t) converge towards zero as t —t 00. 
(Hi) The vector Xi — Xj tends to a limit vector xij, for all i^j = 1, . . . , N. 

In other words, the velocity support collapses exponentially to a single point and the flock 
holds the same disposition. From this theorem we recover the notion of unconditional flocking 
in the regime 7 < \. If 7 > \ in general unconditional flocking doesn't follow, but under some 
conditions on initial data flocking condition is reached, see [13J. 

Note that standard Cucker-Smale model prescribes perfectly symmetric interactions and takes 
in account only the alignment dynamic. As a result total momentum is preserved by the dynamics. 
The introduction of a limited space perception (like a visual cone) breaks symmetry and momentum 
conservation. This choice corresponds to take a function for the strength of the interaction of the 
type 



H-ol V^i? %j 1 ^i ) -tl \\%i %j \)yol y^ii %j 1^1)1 



(2.4) 



3 1 tJ I1 iri ,Xj 3 \) ycx y^i") ^j 1 

where the parameter vector a is related, for example, to the width of the visual cone. 

2.2. D'Orsogna-Bertozzi et al. model. The microscopic model introduced by D'Orsogna, 
Bertozzi et al. [21] considers a self-propelling, attraction and repulsion dynamic and reads 



Vi = (a-b\vi\ 2 )vi - -Tr^2^xiU(\xj -Xi\) 



1,...,A/-, 



(2.5) 



i#i 
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Figure 2.2. The Morse potential in D'Orsogna-Bertozzi et al. model 



where a, b are nonnegative parameters, U : M d — > R is a given potential modeling the short-range 
repulsion and long-range attraction, and N is the number of individuals. Function U gives us the 
attraction-repulsion dynamic typically described by a Morse potential 

r / l A _|_ Ci^p-r/lR 



U(r) = -C A e- r ' lA + C R e~ 



(2.6) 



where Ca, Cr, I a, Ir are positive co nstan ts measuring the strengths and the characteristic lengths 
of the attraction and repulsion. In (2.5) the term (a — b\vi\ 2 )vi characterizes self -propulsion and 
friction. Asymptotically this term give us a desired velocity, in fact for large times the velocity of 
every single particle tends to yjajb. 

The most interesting case in biological applications occurs when the constants in the Morse 
potential satisfy the following inequalities C := Cr/Ca > 1 and I := Ir/Ia < 1? which correspond 
to the long range attraction and short range repulsion. Moreover the choice of the parameters fixes 
the evolution of the N particles system towards a particular equilibrium. The following distinction 
holds: if Cl d > 1 then crystalline patterns are observed and for Cl d < 1 the motion of particles 
converges to a circular motion of constant speed, where d > 2 is the space dimension. In [22J a 
further study of the parameters can be found. 

2.3. Motsch-Tadmor model. In a recent work [35J the authors propose a modification of 
the classical Cucker-Smale model as follows 



X{ — Vi 



N 



N 



J2K x h x j)( v J ~ v i)i 



* = !,.. .,iV, 



(2.7) 



i=i 



where h is defined by 



n\Xi^ %j) — 



H(\xj -Xj 
H(xi) 



1 N 
H(x i ) = -Y,H(\x i -x k \). 



k = l 



The model differs from the classical one, since the influence between two particles H(\xj — xi\) is 
weighted by the average influence on the single individual i. In this way the function h(xi, Xj) lose 
in general any kind of symmetry property of the original Cucker-Smale dynamic. 
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Figure 2.3. One of the possible configurations in the interaction with a perception cone. Individual j is 
perceived by individual i but not vice versa. 

We emphasize, however, that in our general setting this model is included in the Cucker- 



Smale alignment dynamic of type (2.4) with a particular choice for the function defining the space 
perception of the form 

if) a (xi,Xj,Vi) = - . (2.8) 

H(xi) 

This can be interpreted as a higher perception level of zones where the individuals have a higher 
concentration and a lower interest in zones where individuals are more scattered. 

2.4. Perception cone, topological interactions and roosting force. For interacting 
animals like birds, fishes, insects the visual perception of the single individual plays a fundamental 
role [15, 23, 24J. In [13J the authors introduce in the dynamic a further rule: the visual cone. A 
visual cone identifies the area in which interaction is possible and blind area where can not be 
interaction. Mathematically speaking the visual cone depends on an angle, 0, that give us the 
visual width. Together with position and velocity the visual area can be described as follows 

£(^,M) = \yeR d : |fc^4nn ^ cos(0/2)l . (2.9) 

I \{ x i-y)\\ v i\ J 

As already discussed the introduction of a visual cone breaks the typical symmetry of the interaction 



(see Figure 2.3). 



The drawback of this choice is that a single individual that has no one in his visual cone, 
never changes his direction. For real situations this assumption is clearly too strong, since many 
other stimuli are received by the surrounding. We cannot ignore other perceptions like hearing, 
smell and visual memory. For example fishes use their visual perception mostly on large/medium 
distance whereas on medium/short distance they rely on their lateral line. These observations lead 
naturally to improve the idea of a visual cone by introducing a perception cone as follows: we 
assign two different weights measuring the strength/probability of the interaction. A weight p\ in 
the case of strong perception and p2 in case of weak perception, with < p2 < Pi < 1. Note that 
taking pi = 1 and p2 = we have the standard visual cone. For example, in the simulation section 
we consider a perception cone ip a , a = (#,Pi,P2)? with the following form 

l/j a (Xi,Xj,Vi) =p 2 + (pi -P2)^{x ii v ii 0){x j ), (2.10) 
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where ^( Xi ,Vi,e)(') is the indicator function of the set £(o^,i^,0) defined according to (2.9). 

Related efforts to improve the dynamic consider also different ingredients like topological in- 
teractions where individuals interact only with the closest individuals and with a limited number 
of them, see [5j [16]. Another variant concerns the introduction of a term describing a roosting 
force [HJ[T]. In fact, flocking phenomena tends to stay localized in a particular area, this force 
acts orthogonal to the single velocity, giving each particle a tendency towards the origin. 

3. Kinetic equations. For a realistic numerical simulation of a flock the number of inter- 
acting individuals can be rather large, thus we need to solve a very large system of ODEs, which 
can constitute a serious difficulty. An alternative way to tackle this problem is to consider a 
nonnegative distribution function f{x,v,t) describing the number density of individuals at time 
t > in position x e R d with velocity v G M d . The evolution of f(x,v,i) is characterized by a 
kinetic equation which takes into account the motion of individuals due to their own velocity and 
the velocity changes due to the interactions with other individuals. Following [13J we consider 
here binary interaction Boltzmann-type and mean-field kinetic approximation of the microscopic 
dynamics. 



3.1. Boltzman-Povzner kinetic approximations. In agreement with (2.2) and (2.4), we 
consider a microscopic binary interaction between two individuals with positions and velocities 
(x, v) and (y, w) according to 

v* = (1 - rj(H a (\x - y\ , v))v + rjH a (\x - y\ , v)w, 

(3.1) 
w* = rjH a (\x - y\ , w)v + (1 - rj(H a (\x - y\ , w))w, 

where v*,w* are the post-interaction velocities and n a parameter that measures the strength of 
the interaction. Analogous binary interactions can be introduced for other swarming and flocking 
dynamics like D'Orsogna-Bertozzi. 

We describe the interaction of the sistem with following integro- differential equation of Boltz- 
mann type 

(dtf + v ■ Vs/X*, v, t) = -Q(/, f)(x, v, t\ 

r 1 £ (3.2) 

Q(/> /) = / ("f /(^ v ^^)f(v^ ™*,t)~ f(x, v, i)f(y, w, t))dwdy, 

JR 2d J 



where (y*, w*) are the pre- interacting velocity that generate the couple (v, w) according to (3.1 ), J 
is the Jacobian of the transformation of (v, w) to (v*, w*). Without visual limitation the Jacobian 
reads J = (1 — 2r]H(\x — y\)) d . Note that, at variance with classical Boltzmann equation the 
interaction is non local as in Povzner kinetic model [40J. 
Let us introduce the time scaling 

t -> t/e, n = Xe, (3.3) 

where A is a constant and e a small parameter. The scaling corresponds to assume that the 
parameter n characterizing the strength of the microscopic interactions is small, thus the frequency 
of interactions has to increase otherwise the collisional integral will vanish. This corresponds to 
large scale interaction frequencies and small interaction strengths, in agreement with a classical 
mean-field limit and similarly to the so-called grazing collision limit of the Boltzmann equation for 
granular gases [34] . 
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3.2. Derivation of the mean- field kinetic model. First of all let us remark that the 



dynamic (3.1) doesn't preserve the momentum, as consequence of the velocity dependent function 
H a we have 



v* + w* = v + w — rj(H a (\x — y\ ,w) — H a (\x — y\ , v)). 



(3.4) 



Moreover under the assumptions \H a (r,v)\ < 1 and r] < 1/2 , it is easy to prove that the support 
of velocity is limited by initial condition 



v* = (1 -r]H a (\x-y\ ,v))v + rjH a (\x-y\ ,v)w < max{H, \w\}. 



(3.5) 



Considering now the weak formulation of (3.2) the Jacobian term disappears and we get the 
rescaled equation 



d_ 

dt ./ic 



/ (j>(x,v)f(x,Vit)dvdx+ I (v ■ V0(x, v))f(x, v, t)dvdx = 

jR 2d ' ' JR 2d 

x, v*) — <j)(x, v))f(x, v, t)f(y, w, t)dvdxdwdy, 



for t > and for all <\> e C§°(™ 2d 



1 

8 JwL* d 

), such that 



(3.6) 



lim 



4>(x, v)f(x, v, t)dvdx 



R 2d 



4>{x, v)fo(x, v,t)dvdx. 



(3.7) 



where fo(x,v) is the starting density. 

For small values of e we have v* « v thus we can consider the Taylor expansion of (j)(x,v*) 
around v up to the second order we obtain the following formulation to the collisional integral 



(4>(x, v*) — <j)(x, v))f(x, v, t)f(y, w, t)dvdxdwdy = 

= A / (X7 v (f)(x : v) • (w — v))H a (x, y, v)f(x, v, t)f(y, w, t)dvdxdwdy 

JR 4d 



+ A 2 e 



^ d 2 (j)(x,v) 2 



; ,j=l 



dv 2 



(H a (x, y, v)) 2 f(x, v, t)f(y, w, t)dvdxdwdy (3.8) 



■■=h(fj) 



for some v = rv + (1 — r)v*, < r < 1. In the limit e — >• the term ^(Z, f) vanishes since the 
second momentum of the solution is not increasing and H a (x,y,v) < 1 hence [12] 



\h(f,f)\ < 2\\<j,(x,v)\\ CS / \v\ 2 f (x,v)dxdv. 



(3.9) 



Thus in the limit the second-order term can be neglected and ii(/, /) constitutes an approximation 
of the collisional integral Q(/, /), in the strong divergence form 



J r i(/,/) = -v 



v I ( 

JR 2d 



w - v)H a (x, y, v)f(y, w, t)f(x, v, t)dwdy, 



or equivalently in convolution form [12] 

h(f,f) = V„ • {f(x,v,t)[(H a (x,y,v)V v e(v)) * f](x,v,t)} , 



(3.10) 



(3.H) 
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where e(v) = \v\ 2 /2 and * is the (x, ^-convolution. As observed in [12J, the operator h(f,f) 
preserves the dissipation proprieties of original Boltzmann operator. 
Finally we get the mean-field kinetic equation 

dtf + v-V x f = -\V v [£(f)f] (3.12) 

f (/) = / H a (x,y,v)(w-v)f(y,w,t)dwdy. 



As noted in [lj, the continuos kinetic model (3.12) and the microscopic one (2.2)-(2.4) are really 



the same when we take the discrete TV-particle distribution 

N 

N 



1 N 
f(x, v,t) = — ^ S ( x - Xi(t))S(v - Vi(t)), 



where S(-) denotes the Dirac-delta function. 
Remark 3.1. 

• Kinetic formulation for the D'Orsogna Bertozzi et al. with perception cone can be derived 
in the same way and yields the mean-field model 

dtf + v V x f + V v (S(v)f) = -\\7 V • [ / Z a (x, y, v)f(y, w, t)dydw]f(x, v, t), (3.13) 

JR 2d 

where Z a (x,y,v) = [A(x,y) + R(x,y)]i/; a (x,y,v) represents the attraction repulsion term. 

• In 1^3] the authors observed that a certain degree of randomness helps the coherence in the 



collective swarm behavior. Following JTffl . if we add in (2.2) a nonlinear noise term de- 
pending on function H a and perform essentially the same derivation of the above paragraph 
we obtain the kinetic equation 

dtf + v- V x f = -AV V [£(/)/] + aA v (H a * p)f, (3.14) 

where p = p(x, t) represent the mass of the system and a > the strength of the noise. If 
H a (x : y : v) = H(x,y), the right hand side can be written as a Fokker-Plank operator 

V„-(a(ff*p)V„/-A£(/)/), 



and thus a global Maxwellian function is a steady state solution for the equation (3.14) 



3.3. Alternative formulations. In this section we present some alternative formulations 
of the Boltzmann equation describing the binary interaction dynamics for alignment. All the 
formulations share the property that in the mean-field limit originate the same kinetic model 



(3.12). 



The Boltzmann equation (3.2) has much in common with a classical Boltzmann equation for 
Maxwell molecules, in the sense that the collision frequency is independent of the velocity and 
position of individuals. An alternative Boltzmann-like kinetic approximation is obtained with the 
interaction operator 

QUJ)= / H a (x,y,v) \^f(x,v*)f(y,w*)- f(x,v)f(y,w))dwdy, (3.15) 

JR2d \j J 

where now 

V* = (1 — f])v + T]W, 

(3.16) 
w* — rjv + (1 — rj)w. 
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From the modeling viewpoint here the function H a is interpreted as the frequency of interactions 
instead of the strength of the same interactions. 



Clearly the two formulations (3.2) and (3.15) are not equivalent in general. It is easy to verify 



that formally we obtain the same mean- field limit (3.12). Note however that now the second order 



term in the expansion (3.8) is slightly different and reads 

d 2 (j)(x,v) 



Jr 4 






dv 2 



-(Wj-Vj) 



H a (x, y, v)f(x, v, t)f(y, w, t)dvdxdwdy. 



Since H a > (H a ) 2 , in practice we may expect a slower convergence to the mean- field dynamic for 
small values of e. 

Let us finally introduce some stochastic effect in the visual cone perception by defining 

H a (xi,Xj,Vi) = (H(xi,Xj), 

where £ is a random variable distributed accordingly to some b a ((,Xi,Xj,Vi) > s.t. 



®a vS> ' *' 3 '' i) t> ' %iiXjjVi. 



Then the collision term in the form (|3.2|) becomes 

Q(fJ) 



R 2d+1 



6 a (C, x, y, v) [ -jf(x, v*)f(y, w*) - f(x, v)f(y, w) ) dwdyd( 



whereas in the space dependent interaction frequency form (3.15) reads 



Q(f, f)= B a((i x i V> v ) ( -rf( x i v *)f(y, w *) ~ ffa v)f(y, w) 

jR2d+l \J 



dwdydC,, 



(3.17) 



(3.18) 



(3.19) 



where B a (^x,y,v) = b a (£,x,y,v)H(x,y). Again it can be shown that thanks to (3.17) the limit 
asymptotic behavior e — )> is unchanged. We omit the details. 

We conclude the section reporting an example of distribution for the random variable ( which 



corresponds to the stochastic analogue of ( |2.10| ) 

c-- 



with probability pi, 

with probability 1 — pi , 

with probability p2, 

with probability 1 — p2 , 



for y G £(#,?;, 0), 
for ?/ G S(x, v,0), 

for 7/GlR d \S(x,v,6>), 
for 7/GR d \S(x,v,6>). 



4. Monte Carlo methods. Following [21] we introduce different numerical approaches 
for the above kinetic equations based on Monte Carlo methods. The main idea is to approximate 
the dynamic by solving the Boltzmann-like models for small value of e. We will also develop some 



direct Monte Carlo techniques for the limiting kinetic equation (3.12). For the sake of simplicity we 



describe the algorithms in the case of the collision operator (3.2), extensions to the other possible 



formulations presented in Section 3 are also discussed along the section. As we will see, thanks to 
the structure of the equations, the resulting algorithms are fully meshless. 

4.1. Asymptotic binary interaction algorithms. As in most Monte Carlo methods for 
kinetic equations, see [39J, the starting point is a splitting method based on evaluating in two 
different steps the transport and collisional part of the scaled Boltzmann-Povzner equation 



0/ 
dt 



= -v ■ VJ 



(T) 
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% = -Q^f) (C) 



where we used the notation Q e (f, f) to denote the scaled Boltzmann operator (3.2). We emphasize 
that the solution to the collision step for small values of s has very little in common with the classical 
fluid-limit of the Boltzmann equation. Here in fact the whole collision process depends on space 
and on the small scaling parameter e. In particular, in the small e limit the solution is expected 



to converge towards the solution of the mean- field model (3.12). 

By decomposing the collisional operator in equation \C\ in its gain and loss parts we can 
rewrite the collision step as 

% = l[Qi(fj)-pf], (4.i) 

where p > represent the total mass and Qf the gain part of the collisional operator. Without 
loss of generality in the sequel we assume that 

f{x 1 v 1 t)dxdv = 1. 

In order to solve the trasport step we use the exact free flow of the sample particles (xi(t),Vi(t)) 
in a time interval At 

Xi(t + At) = Xi(t) + Vi(t)At, (4.2) 



and thus describe the different schemes used for the interaction process in the form (4.1). 



4.1.1. A Nanbu-like asymptotic method. Let us now consider a time interval [0,T] dis- 
cretized in n to t intervals of size At. We denote by f n the approximation of /(x, v, nAt). 
Thus the forward Euler scheme writes 

/n+l = ^l _ ^j f n + ^Q+(/», /«), (4.3) 

where since f n is a probability density, thanks to mass conservation, also Qf(f n , f n ) is a probability 
density. Under the restriction At < e then also / n+1 is a probability density, since it is a convex 
combination of probability densities. 

From a Monte Carlo point of view equation ( |4.3[ ) can be interpreted as follows: an individual 
with velocity v at position x will not interact with other individuals with probability 1 — At/e 
and it will interact with others with probability At/e according to the interaction law stated by 
Qt(f n i f n )- Since we aim at small values of e the natural choice as in |8J is to take At = e. The 
major difference compare to standard Nanbu algorithm here is the way particles are sampled from 
Qt(f n ->f n ) which does not require the introduction of a space grid. A simple algorithm for the 
solution of ( |4.3[ ) in a time interval [0,T], T = n to t^-, At = s is sketched in the sequel. 

Algorithm 4.1 (Asymptotic Nanbu I). 

1. Given N samples (x^, v®), with k = 1, . . . , N from the initial distribution fo(x, v); 

2. for n = to n to t — 1 

(a) for i = l to N; 

(b) select an index j uniformly among all possible individuals (x^,v^) except i; 

(c) evaluate H a (\xf — x^|, vf); 



(d) compute the velocity change v* using the first relation in (3.1) with r\ = e; 

(e) set(x? +1 ,v? +1 ) = (x?,v*)- 
end for 
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end for 



Next we show how the method extends to the case of collision operator of the type (3.15). In 
this case an acceptance-rejection strategy is used to select interacting individuals since the forward 
Euler scheme reads 

r +i = (i_^ r + ^P+(/», /«), (4.4) 

where P+(/ n , f n ) — Q £ (/,/) + /> is again a probability density. 

Now using the fact that H a < 1 we can adapt the classical acceptance-rejection technique [39] 



to get the following method for (4.4) with At = e 



Algorithm 4.2 (Asymptotic Nanbu II). 

In Algorithm \4.1\ make the following change 
(d) if H a (\xf — x'jIjV^') >£,,£, uniform in [0, 1] then compute the velocity change v* for 
each individual i of pair (i,j) using the first relation in (3.16J) with r] = e; 



(e) set (x'i +1 ,v™ +1 ) = (x™,v*) if the individual has changed its velocity, otherwise set 

Note that in this version two individuals interact always with the same strength in the velocity 
change but with a different probability related to their distance. As a result the total number of 
interactions depend on the distribution of individuals and on average is equal to H a N < N where 

1 N 

Thus the method computes less interactions then the one described in Algorithm |4.1| In fact, 
in regions where individuals are scattered very few interactions will be effectively computed by 
the method. The efficiency of the method can be further improved if one is able to find an easy 
invertible function 1 > W a (xi, Xj,vi) > H a (xi,Xj,Vi) or is capable to compute directly the inverse 
of H a (xi,Xj,Vi). We refer to [39J for further details on these sampling techniques. 

A symmetric version of the previous algorithms which preserves at a microscopic level other 
interaction invariants, like momentum in standard Cucker-Smale model, is obtained as follows 

Algorithm 4.3 (Asymptotic symmetric Nanbu). 

1. Given N samples (x^, vfy, with k = 1, . . . , N from the initial distribution fo(x, v); 

2. for n = to n to t — 1 

(a) set N c = Iround(N/2); 

(b) select N c random pairs (i,j) uniformly without repetition among all possible pairs of 
individuals at time level n. 

(c) evaluate H a (\xf_^-_x]\,vf) and H a (\xf - x^v 7 }); 



relations (3.1) with rj = e; 



(d) For Algorithm J^.l: compute the velocity changes v* , v* for each pair (i,j) using 



(d) For Algorithm 4-2: 



i. 



if H a (\x™ —x^\,v^) > & & uniform in [0, 1] then compute the velocity change v* 



for each pair (i,j) using the first relation in (3.1(fy with r] = e; 



ii. if H a (\xf —x'jljVj') > £j £j uniform in [0, 1] then compute the velocity change v* 



for each pair (i,j) using the second relation in (3.16) with n = e; 



(e) set « +1 ,< +1 ) = (x?,vf), (^ +1 ,^ +1 ) - (£•>*) for all the individuals that 



changed their velocity, 
(f) (x^ +1 ,^ +1 ) = (xfriVft) for all the remaining individuals. 
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end for 
The function Iround(-) denotes the integer stochastic rounding defined as 



Iround(x) 




£<x-[x] 
elsewhere 



where £ is a uniform [0, 1] random number and [•] is the integer part. 

4.1.2. A Bird-like asymptotic method. The most popular Monte Carlo approach to solve 
the collision step in Boltzmann-like equations is due to Bird |7J. The major differences are that the 
method simulate the time continuous equation and that individuals are allowed to interact more 
then once in a single time step. As a result the method achieves a higher time accuracy [39J. 

Here we describe the algorithm for the collision operator described by (3.2). The method is 



based on the observation that the interaction time is a random variable exponentially distributed. 
Thus for N individuals one introduces a local random time counter given by 



At c (0 



N ' 



(4.5) 



with £ a random variable uniformly distributed in [0, 1]. 

A simpler version of the method is based on a constant time counter At c corresponding to the 
average time between interactions. In fact, in a time interval [0, T] we have 



At, 



T 



e 
A' 



(4.6) 



since N c = NT/e is the number of average interactions in the time interval. Of course taking time 
averages the two formulations (4.5) and (4.6) are equivalent. 



From the above considerations, using the symmetric formulation and the time counter At c 
2e/N, we obtain the following method in a time interval [0,T], T = N tot At c 



,N from the initial distribution fo(x,v) 



Algorithm 4.4 (Asymptotic Bird I). 

1. Given N samples (xk,Vk), with k = 1, 

2. forn = to N tot - 1 

(a) select a random pair (i,j) uniformly among all possible pairs; 

(b) evaluate H a (\xi — Xj | , vi) and H a (\x{ — Xj | , Vj ) ; 



(c) compute the velocity changes v* , v* for each pair (i, j) using relations (3.1) with 

V — e 'j 

(d) set Vi - 
end for 



and Vj 



Note that in the above formulation the method has much in common with Algorithm 4.3 except 
for the fact that multiple interactions are allowed during the dynamic (no need to tag particles 
with respect to time level) and that the local time stepping is related to the number of individuals. 
As a result in the limit of large numbers of individuals the method converges towards the time 



continuous Boltzmann equation (3.2) and not to its time discrete counterpart (4.3), as it happens 
for Nanbu formulation. Since in Algorithm 4.3 we have n to t — N tot /N c , the computational cost of 
the methods is the same. 



Similarly Bird's approach can be extended to collision operator in the form (3.15) by intro- 
ducing the following changes 

Algorithm 4.5 (Asymptotic Bird II). 

In Algorithm^. J\ make the following change 
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(c) * if Ha^xf — x'jl^vf) > & £i uniform in [0, 1] then compute the velocity change v\ 
for each pair (i,j) using the first relation in {3. 16) with rj 



* if H a (\xf —x^l^y™) > £j £j uniform in [0, 1] then compute the velocity change v* 



for each pair (i,j) using the second relation in (3.16) with rj = e; 



Finally we sketch the algorithm to implement the stochastic perception cone present in (3.18) 



and (3.19), that can be easily introduced in all the previous algorithms. 



Algorithm 4.6 (Interaction with stochastic perception cone). 

• if Xj E X(#i, Vi, 9) 

— with probability p\ perform the interaction between i and j and compute v[ 
else 

— with probability p2 perform the interaction between i and j and compute v[ 

• if Xi e E(xj,Vj,0) 

— with probability p\ perform the interaction between i and j and compute v'a 
else 

— with probability P2 perform the interaction between i and j and compute v f j 

Note that this reduces further the total number of interactions in the algorithms just described. 
In contrast, for the deterministic case we simply change the relative interaction strengths using 
respectively n = p\S and r] = P2S in the binary interaction rules. 

4.2. Mean- field interaction algorithms. Let us finally tackle directly the limiting mean 
field equation. The interaction step now corresponds to solve 



dtf = -V, 



JR 2d 



f / H a (x, y, v)(v - w)f(y, w, t)dwdy 



As already observed, in a particle setting this corresponds to compute the original 0(N 2 ) dynamic. 
We can reduce the computational cost using a Monte Carlo evaluation of the summation term as 
described in the following simple algorithm. 

Algorithm 4.7 (Mean Field Monte Carlo). 

1. Given N samples v^, with k = 1, . . . , N computed from the initial distribution f$(x, v) and 
M<N; 

2. for n = to n to t ~ 1 

(a) for i = 1 to N 

(b) sample M particles j\ , . . . , Jm uniformly without repetition among all particles; 

(c) compute 

H n (x)- — V# (x- x- v n ) v n - iy ^'^'^V 
oW jlfL a ^' J*'^' * _ m 2-^ H n (x) Jk1 

(d) compute the velocity change 

< +1 = <(1 - AtH2(x t j) + AtH2( Xi )v?. 

end for 
end for 

The overall cost of the above simple algorithm is O = (MN), clearly for M = N we obtain 
the explicit Euler scheme for the original N particle system. In this formulation the method is 
closely related to asymptotic Nanbu's Algorithm |4.1| It is easy to verify that taking M = 1 leads 
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exactly to the same numerical method. On the other hand for M > 1 the above algorithm can be 
interpreted as an averaged asymptotic Nanbu method over M runs since we can rewrite point (d) 
as 



M 



i = l, 



,N. 



fe=i 



The only difference is that averaging the result of Algorithm 4.1 does not guarantee the absence 
of repetitions in the choice of the indexes ji, . . . , Jm- Thus the choice At = e in Algorithm 4.1 
originates a numerical method consistent with the limiting mean-field kinetic equation. Following 
this description we can construct other Monte Carlo methods for the mean field limit taking suitable 
averaged versions of the corresponding algorithms for the Boltzmann models. Here we omit for 
brevity the details. 
Remark 4.1. 

• In Algorithm ^.7| the size of At can be taken larger then the corresponding At = e in 
Algorithm \4.1 However, as we just discussed, since large values of At in the mean-field 
algorithm are essentially equivalent to large values of e in the Boltzmann algorithms we 
don't expect any computational advantage by choosing larger values of At in Algorithm ^. 7[ 

• We remark that changing the time discretization method from Explicit Euler in (4-3) and 



(4.4) to other methods, like semi-implicit methods or method designed for the fluid-limit, 
permits to avoid the stability restriction At < e. Even this approach however does not 
lead to any computational improvement since a strong deterioration in the accuracy of the 
solution is observed when At > e. Here we don't explore further this direction. 

5. Numerical Tests. In this section we first compare accuracy and computational cost of 
some of the different methods and then illustrate their performance on different two-dimensional 
and three-dimensional numerical examples. We use the following notations: ANMC (Algorithm 



4.3), ABMC (Algorithm 4.4), and MFMC M (Algorithm 4.7 for a given M). 



5.1. Accuracy considerations. Here we compare the accuracy of the different algorithms 
studied for a simple space homogeneous situation. In fact, since the algorithms differ only in the 
binary interaction dynamic the homogeneous step is the natural setting for comparing the various 
approaches in term of accuracy. 

We consider the standard Cucker-Smale dynamic. Since we do not have any space dependence 
we assume H(\xi — Xj\) = 1 for each z, j. Thus there is no difference in this test case between 



formulations (3.2) and (3.15) and the relative simulation schemes. 

We take N 
of two gaussian 



50000 individuals and at the initial time the velocity is distributed as the sum 



fo(v) 



(v + vo) 2 
2a?, 



(v ~ vo) 2 \ 



lol 



2tto~ v 



J 



with v = 0.7, a v = a/02. 

The results obtained for ANMC and ABMC with e = 1, 0.1, 0.01 at time T = 1 are reported 
in Figure [5TJ The reference solution is computed using the microscopic model which in this simple 
situation can be solved exactly and gives 



v i (t) = v i (0)e- t + (l-e- t )v, 



^X>(°)- 



j=l 
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Figure 5.1. Convergence to the exact solution (continuous line) of the velocity profiles calculated with ANMC 
(left) and ABMC (right) algorithms. From the top to bottom, At = e with e = 1,0.1,0.01. 



As expected convergence towards the exact solution is observed for both methods. In particular 
for e = 0.01 the results are in good agreement with the direct solution of the microscopic model. 



Next in Figure 5.2 we compare the behavior of the MFMCm method with ANMC for the 
same values of the time step. A considerable difference is observed for large values of At and 
both methods are poorly accurate. On the other hand for smaller values of At they both converge 
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Figure 5.2. MFMCm algorithms compared with ANMC method, at different time steps. From the top 
At = 1, At = 0.1 and At = 0.01. On the left column M = 5 on the left M = 50. 



towards the reference solution. 



Finally in Figure [O] we report the L 2 -norm of the error for ANMC, ABMC and MFMC M 
for various M as a function of At = s. We compare the convergence of the schemes with different 
number of particles TV = 1000 and N = 50000. Note that in both cases the convergence rate of 
the schemes is rather close and for e — At < t*. the statistical error dominates the time error so 
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L„ relative error 



L„ relative error 





Figure 5.3. Relative errors in the L2 norm at T = 1 for the different methods as a function of At = e. On 
the left the error is commuted with N = 1000 particles, on the right the same test is performed with N = 50000 
particles. 



that we observe a saturation effect, where t* « 0.1 for N = 1000 and t* « 0.01 for N = 50000. 

5.2. Computational considerations and ID simulations. Next we want to compare the 
computational cost of the different binary interaction methods for solving the kinetic equation 



(3.12), when compared to the direct numerical solution of the original system (2.1). 



We consider the same one-dimensional test problem as in [13J for the Cucker-Smale dynamic. 
The initial distribution is given by 



fo(x,v) = 



1 



.2(72 



27rcr x a v 



/ -(v + v y 



2(7?, 



\ 



2<72 



for ^o = 2.5, a v = VOA and a x = y/2. 

The computational time for the different methods at T = 1 using e — 0.01 and different 
number of individuals is reported in Table |5.1| Simulations have been performed on a Intel Core 
II dual-core machine using Matlab. The O(N) cost of ANMC and ABMC is evident. The same 
results are also reported in Figure |5.4| which shows the linear growth of the various Monte Carlo 
algorithms. 



N 


10 3 


10 4 


10 5 


10 6 


ANMC 


0.02 s 


0.23 s 


2.82 s 


3.83 x 10 1 s 


ABMC 


0.02 s 


0.21 s 


2.20 s 


3.14 x 10 1 s 


MFMC 50 


0.05 s 


0.41 s 


4.26 s 


4.44 x 10 1 s 


MFMC500 


0.14 s 


1.58 s 


1.33 x 10 1 s 


3.14 x 10 2 s 


MFMC5000 


5.00 s 


5.20 x 10 1 s 


1.71 x 10 3 s 


4.49 x 10 4 s 













Table 5.1 

Computational times for the different methods with various values of N. 
scaling factor e = At is fixed at 10 — 2 . 



The final time is T = 1 and the 
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CPU computation time 




Figure 5.4. Comparison of the computational times for the 
is equal to At = 0.01. 



methods. For each method the time step 



Finally we report the phase-space plots of the previous ID problem obtained using the per- 



ception cone (2.10). Clearly the parameter has no meaning in the one-dimensional case, so that 



the perception limitation concerns only the capability to detect other individuals on the left and 
on the right over the line. 

We compare the evolution in the phase space of two different cases: the classical Cucker-Smale 
model (non visual limitation p 1 = p 2 = 1) and the weighted visual cone with p\ — \ and p 2 = 0.5. 
The results are reported in Figure [575| 

The simulations have been computed using ABMC with At = e = 0.01. The number N of 
individuals is N = 50000, with 7 = 0.05, that is unconditional flock condition. The phase space 
representation is obtained using a space- velocity grid of 100 x 150 cells over the box [—15, 15] x 
[—10, 10]. The results are in good agreement with the one presented in [T3]. Note how perception 
limitations reduce the flocking tendency of the group of individuals by creating two different flocks 
moving towards left and right respectively. 

5.3. 2D Simulations. 

Cucker-Smale dynamics. A generalization of the previous test in two-dimension is obtained 
by considering a group of TV individuals with position (x,y) G M 2 , initially distributed as 



/o(x, y, v x , v y ) = g (x, y)h (v x , v y ), 



where 



9o(x,y) 



27RT 2 



exp{-(x 2 +2/ 2 )/2a 2 }, 



h(r) 



27TV 2 



-(r + ^p) 2 
2z/ 2 



-(r-^o) 2 \ 



2v 2 



with r = \(v x , v y )\, vo = 10, a = y/2 and v = ^/0A. In the following simulations we use N = 100000 
particles and the limited perception cone defined by (2.10). 

We compare the evolution of the space density for different choices of the perception parameters 
and at different times. In the test case considered the parameters for the perception cone are 
= 4/3tt, pi = 1 and p 2 = 0.04. 
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Figure 5.5. ID Cucker-Smale flocking in the phase space. On the left without perception limitations, on the 
right with a perception bound characterized by pi = 1 and p2 = 0.5 



To reconstruct the probability density function in the space we use a 100 x 100 grid on 
[—20,20] x [—20,20]. In each figure we also add the velocity flux to illustrate the flock direc- 
tion. We report the results computed with ABMC method and At = e = 0.01, similar results are 
obtained with the other stochastic binary algorithms. 

At T = 30 the final flocking structure is reached. It is remarkable that in absence of perception 
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Figure 5.6. 2D Cucker-Smale flocking. On the left without perception limitations, on the right with perception 
cone with 9 = 4/37T, p\ = 1 and p2 = 0.04. 



In contrast when we introduce limitations the flocking behavior is less evident and the groups 
splits in two flocks moving in opposite directions. Finally we can also modify the previous example 
to create more complex patterns, but with the same basic structure. 
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Figure 5.7. £D Cucker-Smale dynamics. Spatial density of two flock that merge together. 



The initial distribution now is given by 



g (x, y, v x , v y ) = f (x + m, y, ^, v y ) + / (a; - m, y, v x , v y ), 



where /o is defined as before, and m = 7. We report the results obtained in absence of perception 



cone. The final flocking state is reached at t ~ 30 (see figure 5.7) 



D'Orsogna, Bertozzi model et al. dynamic. Next we want to simulate the D'Orsogna Bertozzi 
model et al. model with the aim to reproduce the typical mill dynamics as in [TTJ [HJ H] but using 
the Boltzmann kinetic approximation. 

Mills and double mills are typical emergence phenomena in school of fishes and flock of birds 
which travel in a compact circular motion, for example, in order to protect themselves from predator 
attacks. At first we work in the twodimensional space taking into account N = 100000 individuals. 



According to the interaction described in (2.5), we consider the long-range attraction and short- 
range repulsion. 



In figure (5.8) the initial data is uniformly distributed on a twodimensional torus, with a 
circular motion. The evolution shows how the attraction and repulsion forces stretch the mill 
reaching after t = 20 a condition of equilibrium in a stable circular motion as a single mill. 
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Figure 5.8. 2D Mill in D'Orsogna-Bertozzi et al. model at various times. Parameters in the attraction- 
repulsion potential are such that C = Cr/Ca = 30, / = Ir/Ia = 0-3, a = 0.7, (3 = 0.05. Final configuration is 
reached after t = 20 



In figure (5.9), we instead consider the following initial data 

(v x +v ) 2 



fo(x,y,v x ,v y ) = 



1 



x 2 + y 2 



47T 2 CT 2 



2a 2 



(V x - ^o) 2 



where a = y/2 and vq = 0.5. Thus density in space is a normal distribution centered in zero and 
velocity distribution has two main directions left and right. The evolution computed with ABMC 
and At = e shows that equilibrium is reached after t = 30 in a stable double mill formation. 

5.4. 3D simulations. Finally we present some three dimensional simulations for the models 
taking into account the different effects of the thee zone dynamic. All the simulations have been 
performed with ABMC and At = e = 0.01. 

Bertozzi-D 'Orsogna et al. model. First we consider the tridimensional extension of the previous 
simulation for the Bertozzi-D 'Orsogna model et al. Initial data is uniformly distributed in space 
on a 3D-torus, and initial velocity is described by a circular motion in the (x, y) components in z 
direction initial velocity has no influence. 
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Figure 5.9. 2D Double Mill in D'Orsogna-Bertozzi et al. model at various times. Parameters in the 
attraction-repulsion potential are such that C = Cr/Ca = 1-6, I = Ir/Ia — 0.025 7 a = 0.7, (3 = 0.05. Final 
configuration is reached after t = 30 



We present the evolution of the swarm mass density and the vectorial field. The equilibrium 
reached after t = 80 is a ball-shaped flock with mass concentrated on the border and empty zones 
in the middle, that is the typical configuration observed for a mill of a fish school. 

Simulation is made taking in account N=200000 particles, and reconstructing the probability 
density function in the space we use a 3D grid with Ax x Ay x Az = 100 x 100 x 100. 



Binary collision algorithms for flocking and swarming 



25 



2. 

0, 








II 




n 



It 



Figure 5.10. Evolution of the 3D mill in D'Orsogna-Bertozzi et al. model at different times. 



Roosting Force. Accordingly to the work [14] we introduce in the D'Orsogna-Bertozzi model 
a roosting force term. The term expresses essentially the tendency of a flock or a school of fishes 
to stay around a certain zone. Such zones usually are of food interest or where birds settle their 
nests. Different approaches can be used to model this biological behavior, see for example [29) 15]. 

Mathematically speaking such term can be described by the introduction of a force term of 
the type 



Froost = ~ [vt • V(/>(>i)] V±. 



(5.1) 



Such force gives the individuals a tendency to move towards the origin, for a suitable function 
Here (/>, called roosting potential is a function <p : R d — > R. In the simulation we take 



<j>(x) 



4\R 



where R r0 ost gives the roosting area radius, and b is a constant weight. Other choice of this roost 
term are of course possible, we refer the interested reader to |14| . 
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Figure 5.11. Trajectory of the center of mass in the roosting dynamic 



Starting from the following initial data 



f (x,y,z,v x ,v y ,v z ) 



4tt 2 



h^ exv {^ [{a 



,/J- 



■ ^o) 2 + (y - yo) 2 + (z- z ) 2 ] + — [vl + v 2 y ] 



2u 2 



with (xo,yo,zo) = (—10, 10,5), after a certain time the simulation shows a flock in stable equilib- 
rium as an orbital motion around the roosting zone. 

The simulation takes in account the following parameters C = Cr/Ca = 30, / = Ir/Ia = 3/5, 
a = 0.7, /3 = 0.05 and the term of roosting force with parameters R r0 ost — 10 and d = 1/10. 
For a long time simulation the center of mass describe the trajectory depicted in figure |5.4| Some 
configurations of the flock at different times obtained using N = 200000 individuals are reported 



in figure 5.12 



6. Conclusion. Mathematical modeling of collective behavior involves the interaction of sev- 
eral individuals (of the order of millions) which may be computationally highly demanding. Here 
we focus on models for flocking and swarming where the particle interactions implies an 0(N 2 ) 
cost for N interacting objects. Using a probabilistic description based on a Boltzmann equation 
we show how it is possible to evaluate the interaction dynamic in only O(N) computations. In 
particular we derive different approximation methods depending on a small parameter e. 

The building block of the method is given by classical binary collision simulations techniques for 
rarefied gas dynamic. Beside the presence of a further scaling parameter the resulting algorithms 
are fully meshless and can be applied to several different microscopic flocking/swarming models. 
Applications of the present ideas to other interacting particle systems and comparison with fast 
multipole methods are under study and will be presented elsewhere. 
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